Twists through turbidity: propagation of light carrying orbital angular momentum through a complex scattering medium

We explore the propagation of structured vortex laser beams-shaped light carrying orbital angular momentum (OAM)-through complex multiple scattering medium. These structured vortex beams consist of a spin component, determined by the polarization of electromagnetic fields, and an orbital component, arising from their spatial structure. Although both spin and orbital angular momenta are conserved when shaped light propagates through a homogeneous, low-scattering medium, we investigate the conservation of these angular momenta during the propagation of Laguerre–Gaussian (LG) beams with varying topological charges through a turbid multiple scattering environment. Our findings demonstrate that the OAM of the LG beam is preserved, exhibiting a distinct phase shift indicative of the ‘twist of light’ through the turbid medium. This preservation of OAM within such environments is confirmed by in-house developed Monte Carlo simulations, showing strong agreement with experimental studies. Our results suggest exciting prospects for leveraging OAM in sensing applications, opening avenues for groundbreaking fundamental research and practical applications in optical communications and remote sensing.


Results
For LG beams, we rely on the following representation valid in paraxial approximation [23][24][25] : Here, k = 2π/ , w(z) = w(0) 1 + (z/z R ) 2 , z R = πw 2 (0)/ , w(0) corresponds to the zero-order Gaussian beam waist and is adjusted to fit experimental image, {ρ, φ, z} is cylindrical coordinate system with beam propagating along z axis.Then phase evolution of the LG beam is defined as 25 : where G(z) = (2p + |ℓ| + 1) arctan(z/z R ) corresponds to the Gouy phase.The helical phase of LG beam, characterized by a phase front that spirals around the beam's axis, results from the azimuthal phase term in the beam's electric field expression.This term correlates with the beam's topological charge ( ℓ and/or radial index (p), defining the beam's unique spatial structure.
The results of computational MC modeling shows that during transmission, the beam not only loses intensity due to phenomena like scattering, refraction, and absorption but also experiences beam spreading (Fig. 1).As one can see the speckle interference pattern, resulting from the superposition of helical wavefront components, exhibits spatial variations both in intensity and phase distributions.Such patterns lead to the distortion of the LG beam's doughnut structure even in low-scattering ( d/l * = 2.5 ) environments (see Fig. 1-(top)) and its complete breakdown in environments characterized by higher scattering ( d/l * = 5).
Here, the complex turbid media, distinguished by either low or multiple light scattering, are defined by their optical depth ( d/l * ), which quantitatively represents the degree of light attenuation and scattering intensity as it travels through the medium; where, d denotes the medium's thickness or the depth of the beam penetration, and l * is the transport mean free path ( l * = 1/µ s , where µ s is the scattering coefficient of the medium).It is generally accepted that d/l * of 5-6 or greater indicates a medium undergoing multiple or diffuse scattering, while values between approximately 2 to 5 suggest single or intermediate scattering, characterized as 'snake-like photons' 26 . (1) LG ℓ p (ρ, φ, z) = 2p!π(|ℓ| + p)!w 2 (z) In contrast, in medium with low scattering ( d/l * ∼ 2 ), the phase of LG beam progresses with minimal dis- tortion, preserving its original OAM state, petal-like helical phase structure (see Fig. 1-(bottom).Conversely, in higher scattering ( d/l * = 5 ), the LG beam's phase structure diffusively spreads, and its helical phase front is disrupted, leading to the formation of a complex phase speckle pattern (refer to Fig. 2).
The numerical OAM sorter employed to confirm presence of the topological charge in the scattered signal is adopted from 27 .It simulates optical elements that, firstly, transform the attenuated helical light into beam with a transverse phase gradient and, secondly, focus the OAM state present in the resulting beam to a specific lateral position defined by the value of topological charge ℓ .This method, originally introduced and experimentally confirmed to efficiently distinguish between OAM states with different charges, exhibits certain potential to detect OAM presence in biomedical diagnostic applications.
Figure 2 demonstrates that the numerical OAM sorter 27 can detect the OAM memory-the retention of the helical phase structure by the LG beam as it propagates through the scattering medium.More specifically we apply a method for efficiently OAM states of light by employing two static optical elements that transform the helical phase structure of OAM beams into a linear phase gradient 27 .This transformation is achieved through a Cartesian to log-polar coordinate mapping, followed by a phase correction to eliminate distortions.The transformed LG beams are then focused to distinct lateral positions corresponding to their OAM states.The horizontal position of the intensity blob(s) on the detector plane (see Fig. 2) is a direct indicator of the OAM state of the detected light.Each unique OAM state results in a distinct and predictable lateral position, allowing for efficient sorting and detection of multiple OAM states simultaneously.
The preservation of the OAM phase, known as OAM phase memory, is ascribed to the congruence in the alterations of the spiral trajectories encircling the beam axis, arising from rotational symmetry.The spiral trajectories undergo substantial modifications from their initial LG beam configurations due to the interaction with the scattering medium.This leads to disparate partial components of the LG beam's helical wavefront undergoing varying extents of phase distortions.Owing to the inherent rotational symmetry, these alterations in the spiral trajectories exhibit uniformity in both magnitude and direction around the beam's axis.Consequently, despite the significant phase distortions resulting from multiple scattering events ( d/l * ∼ 10 ), the overall retention of the LG beam's OAM is distinctly evident in the phase behavior of the speckle pattern.This pattern mirrors the modulation of the LG beam's initial phase as configured by the SLM (see the Supplementary Video).
The results of this phenomenological model are well agreed with the experimental results 28 , presented in Fig. 3.The LG beam preserves its initial phase structure and spatial intensity profile, demonstrating its ability to maintain phase coherence and OAM characteristics in a medium with moderate ( d/l * = 2 ) scattering proper- ties (see Fig. 3a).Propagation through the multiple-scattering medium ( d/l * ∼ 10 ) caused a diffusive spread and disruption of the LG beam's helical phase front, resulting in a complex speckle pattern (see Fig. 3b).Despite significant phase distortions, the phase speckle pattern retained initial modulation, indicating partial OAM memory, especially within the axial annular region (see Fig. 3c,d, respectively).
Thus, despite the prevalence of strong diffuse scattering, the phase within the speckle pattern retains modulation reflective of the LG beam's initial phase, illustrating the persistence of OAM's memory effect amidst multiple scattering.

Discussion
This study embarked on an exploration of the propagation dynamics of LG vortex beams carrying OAM through turbid, complex scattering media.The results presented above illustrate the impact of scattering on the phase preservation of an LG beam as it propagates through different scattering media.In the low-scattering medium ( d/l * = 2 ), the LG beam maintains its helical phase structure, with minimal disruption to its initial phase and intensity profile (see Fig. 3a).This demonstrates that the LG beam can retain its OAM characteristics and phase coherence in environments with moderate scattering properties.Conversely, in the multiple-scattering medium ( d/l * = 9.6 ), the LG beam undergoes significant phase distortion, resulting in a complex speckle pattern (see Fig. 3b).Despite these distortions, the speckle pattern retains some modulation of the initial phase, indicating a degree of OAM memory.This partial preservation of the OAM phase is particularly noticeable within the axial annular region, while it diminishes in the central and outer regions beyond this annular zone (see Fig. 3).
Utilizing both experimental methods and in-house developed MC simulations, we investigated the conservation of OAM and the accompanying phase shifts, indicative of the unique 'twist of light' phenomena through these challenging environments.Our findings demonstrate a robust preservation of OAM despite the multiple scattering processes, underpinning the potential of OAM for enhanced sensing capabilities in optical communications and remote sensing.The employment of a numerical OAM sorter further validated the persistence of OAM in scattered signals, affirming the foundational concept of OAM phase memory within the scattering media.The investigation reveals that despite the inherent scattering challenges presented by complex media, OAM's unique properties remain largely intact, exhibiting only minor distortions even in high-scattering environments.This resilience of OAM through turbid media opens new avenues for advanced sensing applications, offering a novel paradigm for optical communication systems that require high penetration depths and robust signal integrity.The successful implementation of MC simulations aligns closely with experimental observations (Fig. 4), providing a powerful tool for predicting light propagation behaviors in similarly complex environments.This study contributes to the fundamental understanding of light-matter interactions within turbid medium, as well as paves the way for the practical application of OAM in diverse fields such as biomedical imaging, atmospheric sensing, and secure communications.Future work will aim to refine the predictive capabilities of our models and explore the potential of multiplexed OAM states for even more complex communication and sensing tasks, thereby leveraging the full spectrum of opportunities afforded by the unique properties of structured light beams.

Methodology Experimental setup
A Mach-Zehnderinterferometer 29,30 , typically used to characterize the relative phase shift between two laser beams, allows to observe the phase structure of the LG beam interfered with the reference plane wave of the same frequency 31 .In this study, the modified Mach-Zehnder-based interferometer 32 is used to examine an evolution of OAM of the LG beams propagated through the medium (see Fig. 5).
In the experiment the coherent laser source (60 mW, OBIS, Coherent, USA) emitting at 633 nm is used as a light source.To clear the optical mode, laser beam was focused into a single-mode optical fiber (Thorlabs, USA) F. The output beam was collimated using beam collimator (Thorlabs, USA) BC to obtain Gaussian beam with 1.6 mm waist diameter.Also, to obtain the horizontal linear polarization the polarizer (Thorlabs, USA) P was been used after collimator.Further, the Gaussian beam has been split in sample and reference beams using polarizing beam splitter (Thorlabs, USA) PBS.The sample beam illuminated the phase only spatial light modulator (PLUTO-2-NIR-011, Holoeye, Germany) SLM operating in reflective regime.To produce LG beam with different moments the corresponding forked diffraction patterns were generated at SLM.The light diffracted from SLM has been directed using set of mirrors M to the lens L3 ( f = 45 mm , Thorlabs, USA).This lens was used to focus the first order diffraction through pinhole (Thorlabs, USA).Further the LG beam was re-collimated using lens L4 ( f = 45 mm , Thorlabs, USA).Finally, sample beam was passed through the custom made turbid tissue-like media of low scattering: d/l * = 2 (thickness of d = 1 mm and µ s = 10 mm −1 ) and multiple-scattering: d/l * ∼ 10 (thickness of d = 8 mm and µ s = 6 mm −1 ).The anisotropy of scattering in both cases was g = 0.8 .The detailed methodology for the preparation of such phantoms is thoroughly described in the study by Sieryi et al., as referenced in 33 .The reference beam was passed through the half-wave plate (Thorlabs, USA) to control the polarization orientation of Gaussian beam.Further, reference beam was expanded by set of lenses L1 ( f = 30 mm Thorlabs, USA) and L2 ( f = 70 mm , Thorlabs, USA) and directed to the beam splitter (Thorlabs, USA) BS where expanded Gaussian beam interfered with LG beam.Finally, the interference pattern was registered using CMOS camera (DCC3240M, 1280 × 1024 , Thorlabs, USA) CAM in combination with objective (10× , Nikon, Japan) O.The proposed setup allowed to obtain interference patterns in both on-axis 34 and off-axis 35 regimes.The corresponding simulated intensity profiles for Gaussian beam, expanded Gaussian beam, LG 5 0 beam and corresponding interference pattern for on-axis regime are shown in "jet" colour palette at Fig. 5. Black arrows shows the polarization direction.The corresponding simulated phases for LG 5 0 beam at the SLM, before the cuvette with sample liquid and after the cuvette with sample liquid are show in "parula" colour palette at Fig. 5.

Simulation
MC method is vastly used as an effective tool for imitation of image transfer through turbid scattering media 36 , including complex structured ones like biological tissues 37 where the technique is considered as 'gold standard' 38,39 .Recently, a number of various MC methods have been developed to simulate propagation of polarized light within turbid tissue-like scattering media [40][41][42][43] .The MC approach, informed by the exact analytic Milne solution and utilizing the iterative solution to the Bethe-Salpeter equation [44][45][46] with Jones vector formalism, effectively tracks the polarization of MC-photons within turbid tissue-like media and simulates coherent backscattering 47,48 .Advantages of the Bethe-Salpeter-based approach involve a direct relation to the analytic Milne solution 49 and intuitive physical interpretation of multiple scattering process via ladder diagrams.The approach has subsequently progressed to include polarized light, extending the computational capabilities of biomedical optics diagnostic toolbox 50 .Fundamental ground for the polarized MC approaches was established by the Vector Radiative Transfer Equation (VRTE) which represents a system of equations for each Stokes parameter and can be rigorously derived from the Maxwell electromagnetic theory [51][52][53] .Recently, it has been shown on the fundamental level that VRTE and Bethe-Salpeter-based approaches are equivalent under certain conditions 54 .Modern implementations of the polarization-resolved MC 50 aim to provide a comprehensive description of polarized light scattering with either Jones or Mueller formalism, depending on the representation of the polarization state 55 .Recent biomedical polarimetry advancements highlight the efficacy of circularly polarized light (CPL) in characterizing tissues with optical anisotropy [56][57][58] .Accurate simulation tools are essential for studying CPL-tissue interactions 59,60 .
The inclusion of shaped light carrying OAM in MC modeling has garnered attention for its potential to enhance the biomedical diagnosis toolkit by offering deeper insights into the structural complexities of turbid tissue-like media and overcoming existing limitations [61][62][63] .In MC modeling of OAM light, initial phase ψ 0 j is determined according to (2).Correspondingly, initial direction is determined utilizing Poynting vector direction ( � p = {p ρ , p φ , p z }) 23 : where ω is the angular frequency of the light.Regarding the MC method, within the context of corresponding Cartesian coordinates, this is represented by a direction vector ( � s = {s x , s y , s z }): As previously above, the polarization tracing framework, rooted in the iterative solution to the Bethe-Salpeter equation and incorporating both Jones and Stokes-Mueller formalisms 50 , effectively models optical phenomena including reflection and refraction for linear, elliptical, and circular polarizations at the medium's surface.Within the Bethe-Salpeter-based MC model 44,45 , a large amount ( N inc > 10 9 ) of MC-photons with pre-defined statistical weight W j , j = [1 . . .N inc ] is launched from the source oriented under θ i angle to the medium surface, propagates through turbid medium and statistics is collected from those N ph ≤ N inc arrived on the detector.Turbid medium is defined by scattering coefficient µ s , absorption coefficient µ a , anisotropy of scattering g and refractive index n 64 .Each MC-photon defined at the OAM light source is characterized by the initial statistical weight W 0 j , Cartesian coordinates (x 0 j , y 0 j , 0) , propagation direction s 0 j , initial polarization state and, most importantly, by the initial phase ψ 0 j .The topological charge of the LG beam determines both s 0 j and ψ 0 j .To track the polarization state along the trajectory of the MC-photon, we introduce a real-valued vector P , representing the direction of the linearly polarized electric field E 45,46,50 .
After launch, all MC-photons undergo surface ( z = 0 ) interaction and are transmitted to the turbid medium layer with account for Snell's law and appropriate Fresnel coefficients influencing MC-photon weights, directions and polarization.In turbid medium ( z > 0 ) each MC-photon trajectory is modeled as a sequence of the elementary simulations containing limited amount of scattering events N scatt .This procedure has been thoroughly (3)

Scientific Reports
| (2024) 14:20662 | https://doi.org/10.1038/s41598-024-70954-xwww.nature.com/scientificreports/covered in previous works 50,64 .At each i'th scattering event, i = [1 . . .N scatt ] , the following computational steps are performed: random path length l i = −lnξ/µ s is computed (in this paper, we assume that µ a ≪ µ s and ξ ∈ (0, 1] is a uniformly distributed random number), MC-photon is moved to the next position � r i = � r i−1 + � s i l i with weight attenuated according to the Beer-Lambert law ( W i = W i−1 e −µ a l i ), and next propagation direction � s i+1 is evaluated via inversion of the Henyey-Greenstein (HG) phase function 65 where θ ′ is the polar scattering angle in the MC-photon reference plane.Here, we have used position vector � r i = (x i , y i , z i ) and unit direction for each scattering event � s i = [s X , s Y , s Z ] i = [sin θ cos ϕ, sin θ sin ϕ, cos θ] i , with θ, ϕ as azimuthal and polar angles that correspond to the global Cartesian coordinates.It should be noted that, basically, any phase function p can be used 47,66 .If analytical inversion of p is not possible, then table lookup method is involved to ensure fast computational speed.At each step we check if MC-photon path crosses the medium boundary and invoke surface refraction-transmission and detection procedures if this is the case.Evolution of each linearly polarized state � P x = P xx , P xy , P xz , � P y = P yx , P yy , P yz can be traced along MC-photon trajectory r i , i = [1 . . .N scatt ] via the following procedure which is obtained from the iterative solution to Bethe-Salpeter equation 50,63 : where Î is the third-rank unit tensor and ⊗ indicates a direct product.
We repeat outlined computational steps for each scattering event until one of the following conditions is met: either W i < 10 −4 (statistical weight becomes negligible as follows from Beer-Lambert law) or the amount of scattering events N scatt becomes larger than 10 3 .These limitations ensure proper trajectory tracing cut-off 63 .We continue launching MC-photons until the certain amount (no less than N ph = 10 7 ) arrives on the detector.Detection procedure consists of the two checks: MC-photon coordinates should lie within the detector area ( −r d ≤ x N ≤ r d , −r d ≤ y N ≤ r d , z N = 0 ), and refracted direction s N should meet the detector numerical aperture (NA) requirements.We would limit those directions by using acos(s N • s d ) < NA , where s d = [sin(−θ d ), 0, cos(−θ d )] is a unit vector collinear to the detector axis.Both here and in the subsequent sections N is considered to be an index of the detection event.
Thus, each detected MC-photon has defined by the following total statistical weight influenced by its path within medium and depolarization: where 0 < N j < N scatt is the index of detection event for j'th MC-photon, l i is the path length between two neighbouring scattering events and Ŵ R is the Rayleigh factor 45,49,63 .Based on our Bathe-Salpeter based polarization tracking approach, we are also able to introduce polarized W N j and depolarized W ⊥ N j weights for each MC-photon: By ensemble averaging within each detector pixel we are able to obtain both intensity and phase values on the detector: Here, we assume that N px < N ph is the amount of photons that arrived at the specific pixel, N j is the phase of the j-th detected photon, I px is the resulting intensity in this pixel and ψ(x, y) is the resulting phase value in this pixel.
The MC method, a foundational technique in photon transport simulation, has historically encountered the challenge of balancing detailed, accurate modeling with significant computational demands.It has been successfully used for the imitation of image transfer through the complex scattering media, taking into account  36 .The MC is a robust and versatile tool for modeling light propagation through scattering media, valid across various scattering regimes, including Rayleigh ( X/ ≪ 1 ), Mie ( X/ ≈ 1 ), and geometrical optics ( X/ ≫ 1 ) regimes; Here X is the smallest dimension of the main scatterer.Renowned for its exceptional precision in representing the intricate interactions between light and biological tissues, the MC method is widely regarded as the 'gold standard' 39 .Despite this, the MC method's reliance on extensive computational resources has been a longstanding issue, with traditional simulations relying heavily on central processing units (CPUs).While the shift from CPU to GPU processing marked a revolutionary decrease in computation time-transforming processes that once spanned hours into mere seconds-this evolution brought with it an increased burden on power consumption and a heightened environmental impact, reflecting a growing concern in our climate-conscious era.
The above mentioned realization of MC implements, for the first time, the energy-effective MC algorithm for complex light propagation utilizing the recently introduced Apple M-family processors 67 .These cuttingedge chips feature a unified memory architecture, which became a revolutionary advancement in our needs of simulating complex light transfer.This development effectively eliminates the previous constraints on the number of photon packets or statistics that need to be stored and processed, dramatically broadening the scope for MC simulations of photon behaviour.Such capability enables the detailed study/statistical analyses of light propagation through turbid media without the computational and environmental costs historically associated with such efforts.In our evaluations, that the efficiency of these processors has been highlighted by their ability to conduct simulations using approximately 100 times less energy than their traditional GPU counterparts and achieve 300 to 10000 times speed increases than those possible with conventional MC simulations.
The M-family chips possess a sophisticated architecture uniquely optimized for parallel processing, a key factor in their superior performance.We utilize the Metal Compute framework developed by Apple, which plays a crucial role in harnessing these capabilities for MC simulations.Metal provides a low-overhead, hardwareaccelerated graphics and compute Application Programming Interface (API) finely tuned for parallel data processing, offering a C-style programming language that is both powerful and efficient.Within this framework, compute shaders act as programmable kernel functions, specially designed to execute MC algorithms on the GPU, enabling the detailed and complex simulation of photon trajectories with unparalleled efficiency and speed.

Fig. 1 .
Fig. 1.Intensity (top) and phase (bottom) distributions of the LG 5 0 beam along propagation from Spatial Light Modulator (SLM) into the complex scattering medium.Insets show, respectively, 2D spatial distributions of intensity (normalized) and phase at the surface of the medium ( d/l * = 0 ) and within the medium at the depth d/l * = 2.5 and d/l * = 5.0.

Fig. 2 .
Fig. 2. Intensity (a) and phase (b) distributions of the LG 3 0 beam after propagation through a complex scattering medium with 1 mm thickness and specified scattering coefficients µ s = 2, 4, 6, 10 mm −1 ; absorption coefficient µ a = 0.01 mm −1 and anisotropy factor g = 0.8 .Row (c) depicts intensity output of the OAM sorter algorithm adopted from 27 .Dotted red line on left figure indicates position of the intensity peak which corresponds to the Gaussian beam without OAM 27 .The scale bar in (a) represents 0.5 mm, which applies to all images.

Fig. 3 .
Fig. 3. Experimentally observed phase distributions for the LG beam propagating through (a) low-scattering ( d/l * = 2 ) and (b) multiple-scattering ( d/l * = 9.6 ) media.The central annular zone (highlighted by contours) corresponds to the LG 3 0 beam as if it were passing through a non-scattering medium.Phase change at the selected area of the LG 3 0 axial annular area for low-scattering (c), multiple-scattering (d) media and a comparison of relative phase change (e) according to the initial phase configuration ( −3π/10 � 3π/10 ) at the SLM (f).